Simulation of Φ-sat-2 bands from Sentinel-2 bands¶
This notebook implements the simulation chain to derive Φ-sat-2 bands starting from Sentinel-2 L1C images. The processing chain is outlined in the Phisat-2 Mission Overview document.
The aim of this challenge is to create an AI-based application that will be deployed on board of the Φ-sat-2 satellite. In order to train your AI-based models, we simulate Φ-sat-2 imagery starting from Sentinel-2 imagery, which are spectrally very similar. This notebook will allow you to simulate Φ-sat-2 images starting from the entire archive of Sentinel-2 imagery available through Sentinel Hub. In addition, we provide a pre-computed dataset of about 500 Φ-sat-2 images. Checkout Section 13 for more information about the core dataset was generated.
This notebook showcases how to download and simulate the data using eo-learn, a Python library specifically designed to deal with Earth Observation data. The data will be stored into EOPatches as numpy arrays and geopandas dataframes to facilitate processing operations. The output imagery is by default saved in the GeoTiff file format.
The simulation workflow can be summarised as follows:
- download Sentinel-2 L1C bands (
B02,B03,B04,B08,B05,B06,B07) and the Scene ClassificationSCLmask. More details about the bands can be found at this link; - create the cloud, cloud shadow and cirrus masks as separate arrays;
- (optional) filter time-frames by data coverage;
- download metadata about solar irradiance and Earth-Sun distance;
- compute radiances from reflectances;
- compute a pan-chromatic as linear combination of the input bands;
- spatial resampling to Φ-sat-2 pixel size;
- simulate L1A/L1B band-to-band misalignment;
- simulate signal degradation due to Signal-to-Noise Ratio (SNR);
- simulate signal degradation due to Module Transfer Function (MTF);
- for L1C, simulate reflectances from radiance values;
- split the AOI into a smaller set of image chips;
- save bands and masks for the image chips to generate an AI-ready dataset.
Content of the notebook is as follows:
- Requirements
- Data Download
- Metadata Download
- Compute Radiances
- Add Pan-chromatic Band
- Spatial Resampling
- Band Misalignment
- Noise Calculation
- Simulate PSF Filtering
- Alternative Noise and PSF filtering calculation
- Compute L1C
- Gridding
- Export to Tiff
- Create Workflow
- Core Dataset
To run the notebook, install the Python libraries listed in the provided requirements.txt file, for instance by running pip install-r requirements.txt.
%load_ext autoreload
%autoreload 2
%matplotlib inline
import datetime
import os
import matplotlib.pyplot as plt
import cv2
import numpy as np
import geopandas as gpd
from eolearn.core import (
EOTask,
EOPatch,
EOWorkflow,
FeatureType,
MapFeatureTask,
RemoveFeatureTask,
linearly_connect_tasks,
EOExecutor,
)
from eolearn.features import SimpleFilterTask
from eolearn.io import SentinelHubInputTask
from eolearn.features.utils import spatially_resize_image as resize_images
from sentinelhub import (
BBox,
DataCollection,
SHConfig,
get_utm_crs,
wgs84_to_utm,
)
from sentinelhub.exceptions import SHDeprecationWarning
from tqdm.auto import tqdm
from phisat2_constants import (
S2_BANDS,
S2_RESOLUTION,
BBOX_SIZE,
PHISAT2_RESOLUTION,
ProcessingLevels,
WORLD_GDF,
)
from phisat2_utils import (
AddPANBandTask,
AddMetadataTask,
CalculateRadianceTask,
CalculateReflectanceTask,
SCLCloudTask,
BandMisalignmentTask,
PhisatCalculationTask,
AlternativePhisatCalculationTask,
CropTask,
GriddingTask,
ExportGridToTiff,
get_extent,
)
/Users/batic/Work/git/eo/clients/esa/ai4eo/challenge-3/phisat-2/code/phisat2_constants.py:37: FutureWarning: The geopandas.dataset module is deprecated and will be removed in GeoPandas 1.0. You can get the original 'naturalearth_lowres' data from https://www.naturalearthdata.com/downloads/110m-cultural-vectors/.
WORLD_GDF = gpd.read_file(gpd.datasets.get_path("naturalearth_lowres"))
# filter out some SHDeprecationWarnings
import warnings
warnings.filterwarnings("ignore", category=SHDeprecationWarning)
0. Requirements¶
In order to download the training data through the provided APIs you would need to set-up an account for Sentinel Hub.
If you don't have one already, create a free trial account, register an OAuth client from your dashboard by clicking on the "Create a new OAuth client" button, set the "Client grant type" to "Client Credentials" and save the secret value. Paste below the client_id and secret_client_id and you are ready to go. For more information on configuration of SH, please visit this link.
If you need to extend your trial account or have any issue using the service, please contact support at info@sentinel-hub.com.
If you are using the Euro Data Cube (EDC) resources
The accounts have been already created and configured for you, and what you need is only to retrieve the SentinelHub credentials from your dashboard and paste them below (from your EDC dashboard click on EDC Sentinel Hub in the My API service subscriptions, then click on SHOW credential in API Access section).
sh_config = SHConfig()
sh_config.sh_client_id = "<your Sentinel Hub OAuth client ID>"
sh_config.sh_client_secret = "<your Sentinel Hub OAuth client secret>"
Define processing level, AOI and time interval¶
In the following cells you have to specify the processing level, the AOI and the time interval for your simulated Φ-sat-2 bands.
ProcessingLevels._member_names_
['L1A', 'L1B', 'L1C']
PROCESSING_LEVEL = ProcessingLevels.L1C
PROCESSING_LEVEL
<ProcessingLevels.L1C: 'L1C'>
Here the bounding box is computed given the swath width of Φ-sat-2. The user can input the centroid coordinates in WGS84 coordinate reference system, and the bounding box will be automatically created.
def get_utm_bbox(lat_centre: float, lon_centre: float):
"""Returns a bounding box of size corresponding to the swath width of Φ-sat-2 given the centroid of the area-of-interest in WGS84"""
east, north = wgs84_to_utm(lon_centre, lat_centre)
east, north = 10 * int(east / 10), 10 * int(north / 10)
crs = get_utm_crs(lon_centre, lat_centre)
return BBox(
bbox=(
(east - BBOX_SIZE // 2, north - BBOX_SIZE // 2),
(east + BBOX_SIZE // 2, north + BBOX_SIZE // 2),
),
crs=crs,
)
lat_centre, lon_centre = 42.348, 13.397 # l'Aquila
bbox = get_utm_bbox(lat_centre, lon_centre)
bbox
BBox(((357890.0, 4679580.0), (378030.0, 4699720.0)), crs=CRS('32633'))
time_interval = ("2021-09-01", "2021-09-10")
maxcc = 0.7 # The maximum allowed cloud coverage in percent
aux_request_args = {"processing": {"upsampling": "BICUBIC"}}
1. Data Download¶
Sentinel-2 bands and the Scene classification mask are downloaded from Sentinel Hub.
scl_download_task = SentinelHubInputTask(
data_collection=DataCollection.SENTINEL2_L2A,
resolution=S2_RESOLUTION,
additional_data=[(FeatureType.MASK, "SCL")],
maxcc=maxcc,
aux_request_args=aux_request_args,
config=sh_config,
cache_folder="./temp_data/",
time_difference=datetime.timedelta(minutes=180),
)
scl_cloud_task = SCLCloudTask(scl_feature=(FeatureType.MASK, "SCL"))
eop = scl_download_task(bbox=bbox, time_interval=time_interval)
eop = scl_cloud_task(eopatch=eop)
eop
EOPatch(
bbox=BBox(((357890.0, 4679580.0), (378030.0, 4699720.0)), crs=CRS('32633'))
timestamps=[datetime.datetime(2021, 9, 2, 10, 8, 50), datetime.datetime(2021, 9, 7, 10, 8, 54)]
mask={
SCL_CIRRUS: numpy.ndarray(shape=(2, 2014, 2014, 1), dtype=uint8)
SCL_CLOUD: numpy.ndarray(shape=(2, 2014, 2014, 1), dtype=uint8)
SCL_CLOUD_SHADOW: numpy.ndarray(shape=(2, 2014, 2014, 1), dtype=uint8)
}
)
additional_bands = [(FeatureType.DATA, name) for name in ["sunZenithAngles"]]
masks = [(FeatureType.MASK, "dataMask")]
input_task = SentinelHubInputTask(
data_collection=DataCollection.SENTINEL2_L1C,
resolution=S2_RESOLUTION,
bands_feature=(FeatureType.DATA, "BANDS"),
additional_data=masks + additional_bands,
bands=S2_BANDS, # note the order of these bands, where B08 follows B03
aux_request_args=aux_request_args,
config=sh_config,
cache_folder="./temp_data/",
time_difference=datetime.timedelta(minutes=180),
)
eop = input_task(eopatch=eop)
eop
EOPatch(
bbox=BBox(((357890.0, 4679580.0), (378030.0, 4699720.0)), crs=CRS('32633'))
timestamps=[datetime.datetime(2021, 9, 2, 10, 8, 50), datetime.datetime(2021, 9, 7, 10, 8, 54)]
data={
BANDS: numpy.ndarray(shape=(2, 2014, 2014, 7), dtype=float32)
sunZenithAngles: numpy.ndarray(shape=(2, 2014, 2014, 1), dtype=float32)
}
mask={
SCL_CIRRUS: numpy.ndarray(shape=(2, 2014, 2014, 1), dtype=uint8)
SCL_CLOUD: numpy.ndarray(shape=(2, 2014, 2014, 1), dtype=uint8)
SCL_CLOUD_SHADOW: numpy.ndarray(shape=(2, 2014, 2014, 1), dtype=uint8)
dataMask: numpy.ndarray(shape=(2, 2014, 2014, 1), dtype=bool)
}
)
Visualise the downloaded bands as well as the cloud, cloud shadow and cirrus masks.
n_timestamps = len(eop.timestamp)
_, axes = plt.subplots(ncols=2, nrows=4, figsize=(10, 20), sharex=True, sharey=True)
for t_idx, timestamp in enumerate(eop.timestamp):
axes[0][t_idx].imshow(np.clip(2.5 * eop.data["BANDS"][t_idx][..., [2, 1, 0]], 0, 1))
axes[0][t_idx].set_title(eop.timestamp[t_idx])
axes[1][t_idx].imshow(eop.mask["SCL_CLOUD"][t_idx][..., 0], vmin=0, vmax=1)
axes[1][t_idx].set_title("SCL_CLOUD")
axes[2][t_idx].imshow(eop.mask["SCL_CIRRUS"][t_idx][..., 0], vmin=0, vmax=1)
axes[2][t_idx].set_title("SCL_CIRRUS")
axes[3][t_idx].imshow(eop.mask["SCL_CLOUD_SHADOW"][t_idx][..., 0], vmin=0, vmax=1)
axes[3][t_idx].set_title("SCL_CLOUD_SHADOW")
plt.tight_layout()
/var/folders/wv/d2zg_qp156z_4hsjy3yb2shm0000gn/T/ipykernel_95471/1616319318.py:1: EODeprecationWarning: The attribute `timestamp` is deprecated, use `timestamps` instead. n_timestamps = len(eop.timestamp) /var/folders/wv/d2zg_qp156z_4hsjy3yb2shm0000gn/T/ipykernel_95471/1616319318.py:4: EODeprecationWarning: The attribute `timestamp` is deprecated, use `timestamps` instead. for t_idx, timestamp in enumerate(eop.timestamp): /var/folders/wv/d2zg_qp156z_4hsjy3yb2shm0000gn/T/ipykernel_95471/1616319318.py:6: EODeprecationWarning: The attribute `timestamp` is deprecated, use `timestamps` instead. axes[0][t_idx].set_title(eop.timestamp[t_idx])
(Optional) Filter by valid data¶
Depending on the location and time interval selected, no data can be available. This task removes the time-frames which have no data.
The same task can be used to filter by cloud coverage. Simply create a counter-part of the full_valid_data method with the desired threshold, and adapt the SimpleFilterTask to use the cloud mask and the newly created method.
def full_valid_data(array: np.array) -> bool:
return np.mean(array) == 1.0
filter_task = SimpleFilterTask((FeatureType.MASK, "dataMask"), full_valid_data)
eop = filter_task(eop)
eop
EOPatch(
bbox=BBox(((357890.0, 4679580.0), (378030.0, 4699720.0)), crs=CRS('32633'))
timestamps=[datetime.datetime(2021, 9, 2, 10, 8, 50), datetime.datetime(2021, 9, 7, 10, 8, 54)]
data={
BANDS: numpy.ndarray(shape=(2, 2014, 2014, 7), dtype=float32)
sunZenithAngles: numpy.ndarray(shape=(2, 2014, 2014, 1), dtype=float32)
}
mask={
SCL_CIRRUS: numpy.ndarray(shape=(2, 2014, 2014, 1), dtype=uint8)
SCL_CLOUD: numpy.ndarray(shape=(2, 2014, 2014, 1), dtype=uint8)
SCL_CLOUD_SHADOW: numpy.ndarray(shape=(2, 2014, 2014, 1), dtype=uint8)
dataMask: numpy.ndarray(shape=(2, 2014, 2014, 1), dtype=bool)
}
)
n_timestamps = len(eop.timestamp)
_, axes = plt.subplots(ncols=2, nrows=1, figsize=(15, 15), sharex=True, sharey=True)
for t_idx, timestamp in enumerate(eop.timestamp):
axes[t_idx].imshow(np.clip(2.5 * eop.data["BANDS"][t_idx][..., [2, 1, 0]], 0, 1))
axes[t_idx].set_title(eop.timestamp[t_idx])
plt.tight_layout()
/var/folders/wv/d2zg_qp156z_4hsjy3yb2shm0000gn/T/ipykernel_95471/95940776.py:1: EODeprecationWarning: The attribute `timestamp` is deprecated, use `timestamps` instead. n_timestamps = len(eop.timestamp) /var/folders/wv/d2zg_qp156z_4hsjy3yb2shm0000gn/T/ipykernel_95471/95940776.py:4: EODeprecationWarning: The attribute `timestamp` is deprecated, use `timestamps` instead. for t_idx, timestamp in enumerate(eop.timestamp): /var/folders/wv/d2zg_qp156z_4hsjy3yb2shm0000gn/T/ipykernel_95471/95940776.py:6: EODeprecationWarning: The attribute `timestamp` is deprecated, use `timestamps` instead. axes[t_idx].set_title(eop.timestamp[t_idx])
2. Metadata Download¶
Download the metadata necessary to compute radiance values. Such metadata consists of the solar irradiance per band, and the Earth-Sun distance.
add_meta_task = AddMetadataTask(config=sh_config)
eop = add_meta_task(eop)
/Users/batic/Work/git/eo/clients/esa/ai4eo/challenge-3/phisat-2/code/phisat2_utils.py:205: EODeprecationWarning: The attribute `timestamp` is deprecated, use `timestamps` instead. if not all([eopatch, eopatch.bbox, eopatch.timestamp]): /Users/batic/Work/git/eo/clients/esa/ai4eo/challenge-3/phisat-2/code/phisat2_utils.py:215: EODeprecationWarning: The attribute `timestamp` is deprecated, use `timestamps` instead. time=[eopatch.timestamp[0], eopatch.timestamp[-1]], /Users/batic/Work/git/eo/clients/esa/ai4eo/challenge-3/phisat-2/code/phisat2_utils.py:217: EODeprecationWarning: The attribute `timestamp` is deprecated, use `timestamps` instead. tiles = self.filter_and_sort_tiles(list(query), eopatch.timestamp) /Users/batic/Work/git/eo/clients/esa/ai4eo/challenge-3/phisat-2/code/phisat2_utils.py:219: EODeprecationWarning: The attribute `timestamp` is deprecated, use `timestamps` instead. dim = len(eopatch.timestamp) /Users/batic/Work/git/eo/clients/esa/ai4eo/challenge-3/phisat-2/code/phisat2_utils.py:224: EODeprecationWarning: The attribute `timestamp` is deprecated, use `timestamps` instead. for tile_idx, (tile, timestamp) in enumerate(zip(tiles, eopatch.timestamp)):
# check result:
eop
EOPatch(
bbox=BBox(((357890.0, 4679580.0), (378030.0, 4699720.0)), crs=CRS('32633'))
timestamps=[datetime.datetime(2021, 9, 2, 10, 8, 50), datetime.datetime(2021, 9, 7, 10, 8, 54)]
data={
BANDS: numpy.ndarray(shape=(2, 2014, 2014, 7), dtype=float32)
sunZenithAngles: numpy.ndarray(shape=(2, 2014, 2014, 1), dtype=float32)
}
mask={
SCL_CIRRUS: numpy.ndarray(shape=(2, 2014, 2014, 1), dtype=uint8)
SCL_CLOUD: numpy.ndarray(shape=(2, 2014, 2014, 1), dtype=uint8)
SCL_CLOUD_SHADOW: numpy.ndarray(shape=(2, 2014, 2014, 1), dtype=uint8)
dataMask: numpy.ndarray(shape=(2, 2014, 2014, 1), dtype=bool)
}
scalar={
earth_sun_dist: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B02: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B03: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B04: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B05: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B06: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B07: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B08: numpy.ndarray(shape=(2, 1), dtype=float64)
}
)
# plot sol_irr values for each band:
fig, ax = plt.subplots(figsize=(12, 6))
for band in S2_BANDS:
ax.plot_date(
x=eop.timestamp,
y=eop.scalar[f"sol_irr_{band}"],
linestyle="solid",
label=f"sol_irr_{band}",
)
fig.legend()
ax.grid()
/var/folders/wv/d2zg_qp156z_4hsjy3yb2shm0000gn/T/ipykernel_95471/1309170822.py:6: EODeprecationWarning: The attribute `timestamp` is deprecated, use `timestamps` instead. x=eop.timestamp,
3. Compute radiances¶
Radiances are computed from the original reflectance values. The following processing steps, e.g. SNR simulation and PSF filtering, are valid if applied to radiances.
radiance_task = CalculateRadianceTask(
(FeatureType.DATA, "BANDS"), (FeatureType.DATA, "BANDS-RAD")
)
eop = radiance_task(eop)
eop
EOPatch(
bbox=BBox(((357890.0, 4679580.0), (378030.0, 4699720.0)), crs=CRS('32633'))
timestamps=[datetime.datetime(2021, 9, 2, 10, 8, 50), datetime.datetime(2021, 9, 7, 10, 8, 54)]
data={
BANDS: numpy.ndarray(shape=(2, 2014, 2014, 7), dtype=float32)
BANDS-RAD: numpy.ndarray(shape=(2, 2014, 2014, 7), dtype=float64)
sunZenithAngles: numpy.ndarray(shape=(2, 2014, 2014, 1), dtype=float32)
}
mask={
SCL_CIRRUS: numpy.ndarray(shape=(2, 2014, 2014, 1), dtype=uint8)
SCL_CLOUD: numpy.ndarray(shape=(2, 2014, 2014, 1), dtype=uint8)
SCL_CLOUD_SHADOW: numpy.ndarray(shape=(2, 2014, 2014, 1), dtype=uint8)
dataMask: numpy.ndarray(shape=(2, 2014, 2014, 1), dtype=bool)
}
scalar={
earth_sun_dist: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B02: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B03: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B04: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B05: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B06: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B07: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B08: numpy.ndarray(shape=(2, 1), dtype=float64)
}
)
t_idx = 0
_, axes = plt.subplots(ncols=2, nrows=1, figsize=(15, 15), sharex=True, sharey=True)
axes[0].imshow(np.clip(2.5 * eop.data["BANDS"][t_idx][..., [2, 1, 0]], 0, 1))
axes[0].set_title(f"REFLECTANCES - {eop.timestamp[t_idx]}")
axes[1].imshow(np.clip(eop.data["BANDS-RAD"][t_idx][..., [2, 1, 0]] / 100, 0, 1))
axes[1].set_title(f"RADIANCES - {eop.timestamp[t_idx]}")
plt.tight_layout()
/var/folders/wv/d2zg_qp156z_4hsjy3yb2shm0000gn/T/ipykernel_95471/811664563.py:5: EODeprecationWarning: The attribute `timestamp` is deprecated, use `timestamps` instead.
axes[0].set_title(f"REFLECTANCES - {eop.timestamp[t_idx]}")
/var/folders/wv/d2zg_qp156z_4hsjy3yb2shm0000gn/T/ipykernel_95471/811664563.py:8: EODeprecationWarning: The attribute `timestamp` is deprecated, use `timestamps` instead.
axes[1].set_title(f"RADIANCES - {eop.timestamp[t_idx]}")
4. Add pan-chromatic band¶
Compute a pan-chromatic band using a weighted linear combination of the input bands. The pan-chromatic band is inserted at index 3, in line with the acquisition order of the Φ-sat-2 payload.
add_pan_task = AddPANBandTask(
(FeatureType.DATA, "BANDS-RAD"), (FeatureType.DATA, "BANDS-RAD-PAN")
)
eop = add_pan_task(eop)
eop
EOPatch(
bbox=BBox(((357890.0, 4679580.0), (378030.0, 4699720.0)), crs=CRS('32633'))
timestamps=[datetime.datetime(2021, 9, 2, 10, 8, 50), datetime.datetime(2021, 9, 7, 10, 8, 54)]
data={
BANDS: numpy.ndarray(shape=(2, 2014, 2014, 7), dtype=float32)
BANDS-RAD: numpy.ndarray(shape=(2, 2014, 2014, 7), dtype=float64)
BANDS-RAD-PAN: numpy.ndarray(shape=(2, 2014, 2014, 8), dtype=float64)
sunZenithAngles: numpy.ndarray(shape=(2, 2014, 2014, 1), dtype=float32)
}
mask={
SCL_CIRRUS: numpy.ndarray(shape=(2, 2014, 2014, 1), dtype=uint8)
SCL_CLOUD: numpy.ndarray(shape=(2, 2014, 2014, 1), dtype=uint8)
SCL_CLOUD_SHADOW: numpy.ndarray(shape=(2, 2014, 2014, 1), dtype=uint8)
dataMask: numpy.ndarray(shape=(2, 2014, 2014, 1), dtype=bool)
}
scalar={
earth_sun_dist: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B02: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B03: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B04: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B05: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B06: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B07: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B08: numpy.ndarray(shape=(2, 1), dtype=float64)
}
)
5. Spatial Resampling¶
Spatially resample the original Sentinel-2 pixel size of 10 m to Φ-sat-2 pixel size of 4.75 m.
features_to_resize = {
FeatureType.DATA: ["BANDS-RAD-PAN", "sunZenithAngles"],
FeatureType.MASK: [
"SCL_CLOUD",
"SCL_CIRRUS",
"SCL_CLOUD_SHADOW",
"dataMask",
],
}
NEW_SIZE = (int(BBOX_SIZE / PHISAT2_RESOLUTION), int(BBOX_SIZE / PHISAT2_RESOLUTION))
# covert dataMask to uint8
casting_task = MapFeatureTask(
(FeatureType.MASK, "dataMask"), (FeatureType.MASK, "dataMask"), np.uint8
)
eop = casting_task(eop)
resize_task_list = []
for feature_type in tqdm(features_to_resize.keys()):
for feature in tqdm(features_to_resize[feature_type]):
resize_task_list.append(
MapFeatureTask(
(feature_type, feature),
(feature_type, f"{feature}_RES"),
resize_images,
new_size=NEW_SIZE,
resize_method="nearest",
)
)
eop = resize_task_list[-1](eop)
0%| | 0/2 [00:00<?, ?it/s]
0%| | 0/2 [00:00<?, ?it/s]
0%| | 0/4 [00:00<?, ?it/s]
Visualise resampled bands. Note the size of the resampled arrays.
t_idx = 1
fig, axs = plt.subplots(figsize=(15, 15))
axs.imshow(eop.data["BANDS-RAD-PAN_RES"][t_idx][..., [2, 1, 0]] / 100);
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Remove unnecessary features¶
This step removes from memory arrays that are no longer needed.
remove_feature_task1 = RemoveFeatureTask(
[
(FeatureType.DATA, "BANDS"),
(FeatureType.DATA, "BANDS-RAD"),
(FeatureType.DATA, "BANDS-RAD-PAN"),
(FeatureType.DATA, "sunZenithAngles"),
(FeatureType.MASK, "SCL_CLOUD"),
(FeatureType.MASK, "SCL_CLOUD_SHADOW"),
(FeatureType.MASK, "SCL_CIRRUS"),
(FeatureType.MASK, "dataMask"),
]
)
eop = remove_feature_task1(eop)
eop
EOPatch(
bbox=BBox(((357890.0, 4679580.0), (378030.0, 4699720.0)), crs=CRS('32633'))
timestamps=[datetime.datetime(2021, 9, 2, 10, 8, 50), datetime.datetime(2021, 9, 7, 10, 8, 54)]
data={
BANDS-RAD-PAN_RES: numpy.ndarray(shape=(2, 4240, 4240, 8), dtype=float64)
sunZenithAngles_RES: numpy.ndarray(shape=(2, 4240, 4240, 1), dtype=float32)
}
mask={
SCL_CIRRUS_RES: numpy.ndarray(shape=(2, 4240, 4240, 1), dtype=uint8)
SCL_CLOUD_RES: numpy.ndarray(shape=(2, 4240, 4240, 1), dtype=uint8)
SCL_CLOUD_SHADOW_RES: numpy.ndarray(shape=(2, 4240, 4240, 1), dtype=uint8)
dataMask_RES: numpy.ndarray(shape=(2, 4240, 4240, 1), dtype=uint8)
}
scalar={
earth_sun_dist: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B02: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B03: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B04: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B05: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B06: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B07: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B08: numpy.ndarray(shape=(2, 1), dtype=float64)
}
)
6. Band misalignment¶
The simulation of three processing levels is supported, namely L1A, L1B and L1C. L1A differs in the band misalignment process compared to L1B and L1C, with L1B and L1C displaying less misalignment compared to L1A.
L1A processing level:
- define a constant with shifts wrt reference Φ-sat-2 B6 band
np.array([0,1.105,1.046,0.943,0.837,0.717,0.55,0.44,]); - draw random shifts from N(0.0, 0.4)
np.random.normal(0.0, 0.4, size=(8,)); - calculate the shift vector magnitudes by adding the random and constant shifts;
- draw random angles fi from (0 - 2Pi) for each band;
- multiply the shifts to the random angle vector for each band;
- accummulate shifts accross bands;
- apply the
cv2.warpAffinemethod looping through the channels, and applying the random shifts computed above- use 0 as filling value at the border;
- use Nearest Neighboord interpolation to avoid any additional smoothing at this stage.
L1B/L1C processing levels:
- draw an amplitude value following a normal distribution of N(0,1) over land and N(0,6) over a full ocean surface scene. If the scene contained both land/sea, N(0,1) is used by default. The possibility to use N(0,6) should be a user-defined parametrization;
- draw a fully random angle from 0 – 2pi;
- compute the associated X and Y from the amplitude/angle;
- apply this procedure independently for all bands except for the Red band, which is considered as the reference band.
band_shift_task = BandMisalignmentTask(
(FeatureType.DATA, "BANDS-RAD-PAN_RES"),
(FeatureType.DATA, "BANDS-RAD-PAN_SHIFTED"),
PROCESSING_LEVEL,
std_sea=6,
interpolation_method=cv2.INTER_NEAREST,
)
eop = band_shift_task(eop)
# check result
eop
EOPatch(
bbox=BBox(((357890.0, 4679580.0), (378030.0, 4699720.0)), crs=CRS('32633'))
timestamps=[datetime.datetime(2021, 9, 2, 10, 8, 50), datetime.datetime(2021, 9, 7, 10, 8, 54)]
meta_info={
Shifts: {0: [(-0.1312835857123827, 0.6266645003616177), ...]<length=8>, ...}<length=2>
}
data={
BANDS-RAD-PAN_RES: numpy.ndarray(shape=(2, 4240, 4240, 8), dtype=float64)
BANDS-RAD-PAN_SHIFTED: numpy.ndarray(shape=(2, 4240, 4240, 8), dtype=float64)
sunZenithAngles_RES: numpy.ndarray(shape=(2, 4240, 4240, 1), dtype=float32)
}
mask={
SCL_CIRRUS_RES: numpy.ndarray(shape=(2, 4240, 4240, 1), dtype=uint8)
SCL_CLOUD_RES: numpy.ndarray(shape=(2, 4240, 4240, 1), dtype=uint8)
SCL_CLOUD_SHADOW_RES: numpy.ndarray(shape=(2, 4240, 4240, 1), dtype=uint8)
dataMask_RES: numpy.ndarray(shape=(2, 4240, 4240, 1), dtype=uint8)
}
scalar={
earth_sun_dist: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B02: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B03: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B04: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B05: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B06: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B07: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B08: numpy.ndarray(shape=(2, 1), dtype=float64)
}
)
# With Nearest Neighbour interpolation, the effective shifts applied are rounded to integer values.
for shifts in eop.meta_info["Shifts"].values():
print(np.round(shifts))
[[-0. 1.] [ 1. 1.] [ 0. 0.] [-0. -1.] [-0. 0.] [-1. -0.] [-1. 1.] [-0. -0.]] [[ 2. 0.] [ 0. -0.] [ 0. 0.] [-0. 0.] [-0. 0.] [-0. -2.] [-1. 1.] [ 1. -0.]]
# plot the shift in bands
ts_idx = 1
num_bands = 8
for b_idx in range(0, num_bands):
_, ax = plt.subplots(ncols=3, nrows=1, figsize=(15, 15), sharex=True, sharey=True)
ax[0].imshow(
eop.data["BANDS-RAD-PAN_RES"][ts_idx][-100:, :100, b_idx] / 100, vmin=0, vmax=1
)
ax[0].set_title(f"Original band - index {b_idx}")
ax[1].imshow(
eop.data["BANDS-RAD-PAN_SHIFTED"][ts_idx][-100:, :100, b_idx] / 100,
vmin=0,
vmax=1,
)
ax[1].set_title(f"Shifted band- index {b_idx}")
diff = (
eop.data["BANDS-RAD-PAN_RES"][ts_idx][-100:, :100, b_idx]
- eop.data["BANDS-RAD-PAN_SHIFTED"][ts_idx][-100:, :100, b_idx]
)
ax[2].imshow(diff)
ax[2].set_title("Difference")
remove_feature_task2 = RemoveFeatureTask([(FeatureType.DATA, "BANDS-RAD-PAN_RES")])
eop = remove_feature_task2(eop)
eop
EOPatch(
bbox=BBox(((357890.0, 4679580.0), (378030.0, 4699720.0)), crs=CRS('32633'))
timestamps=[datetime.datetime(2021, 9, 2, 10, 8, 50), datetime.datetime(2021, 9, 7, 10, 8, 54)]
meta_info={
Shifts: {0: [(-0.1312835857123827, 0.6266645003616177), ...]<length=8>, ...}<length=2>
}
data={
BANDS-RAD-PAN_SHIFTED: numpy.ndarray(shape=(2, 4240, 4240, 8), dtype=float64)
sunZenithAngles_RES: numpy.ndarray(shape=(2, 4240, 4240, 1), dtype=float32)
}
mask={
SCL_CIRRUS_RES: numpy.ndarray(shape=(2, 4240, 4240, 1), dtype=uint8)
SCL_CLOUD_RES: numpy.ndarray(shape=(2, 4240, 4240, 1), dtype=uint8)
SCL_CLOUD_SHADOW_RES: numpy.ndarray(shape=(2, 4240, 4240, 1), dtype=uint8)
dataMask_RES: numpy.ndarray(shape=(2, 4240, 4240, 1), dtype=uint8)
}
scalar={
earth_sun_dist: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B02: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B03: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B04: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B05: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B06: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B07: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B08: numpy.ndarray(shape=(2, 1), dtype=float64)
}
)
Crop the image to contain only valid data after band misalignment
crop_task = CropTask(
features_to_crop=[
(FeatureType.DATA, "BANDS-RAD-PAN_SHIFTED"),
(FeatureType.DATA, "sunZenithAngles_RES"),
(FeatureType.MASK, "SCL_CLOUD_RES"),
(FeatureType.MASK, "SCL_CLOUD_SHADOW_RES"),
(FeatureType.MASK, "SCL_CIRRUS_RES"),
(FeatureType.MASK, "dataMask_RES"),
]
)
eop = crop_task(eop)
eop
EOPatch(
bbox=BBox(((358232.0, 4679922.0), (377688.0, 4699378.0)), crs=CRS('32633'))
timestamps=[datetime.datetime(2021, 9, 2, 10, 8, 50), datetime.datetime(2021, 9, 7, 10, 8, 54)]
meta_info={
Shifts: {0: [(-0.1312835857123827, 0.6266645003616177), ...]<length=8>, ...}<length=2>
}
data={
BANDS-RAD-PAN_SHIFTED: numpy.ndarray(shape=(2, 4096, 4096, 8), dtype=float64)
sunZenithAngles_RES: numpy.ndarray(shape=(2, 4096, 4096, 1), dtype=float32)
}
mask={
SCL_CIRRUS_RES: numpy.ndarray(shape=(2, 4096, 4096, 1), dtype=uint8)
SCL_CLOUD_RES: numpy.ndarray(shape=(2, 4096, 4096, 1), dtype=uint8)
SCL_CLOUD_SHADOW_RES: numpy.ndarray(shape=(2, 4096, 4096, 1), dtype=uint8)
dataMask_RES: numpy.ndarray(shape=(2, 4096, 4096, 1), dtype=uint8)
}
scalar={
earth_sun_dist: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B02: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B03: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B04: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B05: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B06: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B07: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B08: numpy.ndarray(shape=(2, 1), dtype=float64)
}
)
7. Noise Calculation¶
This step simulates noise according to the signal-toNoise (SNR) specifications of the optical system.
A compiled executable will perform the SNR and PSF calculation.
Executables for the following Operating Systems (OSs) have been compiled:
- Unix,
phisat2_unix.bin - Windows,
phisat2_win.bin - MacOS,
phisat2_osx-arm64.binfor ARM chips,phisat2_osx-x86_64.binfor Intel chips
Download the suitable binary for your OS from this link.
Once downloaded, allow the operating system to run the file and make it executable (e.g. chmod a+x phisat2_unix.bin). Then set the path to the executable in the cell below.
EXEC_PATH = "executables/phisat2_osx-arm64.bin"
snr_task = PhisatCalculationTask(
input_feature=(FeatureType.DATA, "BANDS-RAD-PAN_SHIFTED"),
output_feature=(FeatureType.DATA, "L_out_RES"),
executable=EXEC_PATH,
calculation="SNR",
)
eop = snr_task.execute(eop)
# check
eop
EOPatch(
bbox=BBox(((358232.0, 4679922.0), (377688.0, 4699378.0)), crs=CRS('32633'))
timestamps=[datetime.datetime(2021, 9, 2, 10, 8, 50), datetime.datetime(2021, 9, 7, 10, 8, 54)]
meta_info={
Shifts: {0: [(-0.1312835857123827, 0.6266645003616177), ...]<length=8>, ...}<length=2>
}
data={
BANDS-RAD-PAN_SHIFTED: numpy.ndarray(shape=(2, 4096, 4096, 8), dtype=float64)
L_out_RES: numpy.ndarray(shape=(2, 4096, 4096, 8), dtype=float64)
sunZenithAngles_RES: numpy.ndarray(shape=(2, 4096, 4096, 1), dtype=float32)
}
mask={
SCL_CIRRUS_RES: numpy.ndarray(shape=(2, 4096, 4096, 1), dtype=uint8)
SCL_CLOUD_RES: numpy.ndarray(shape=(2, 4096, 4096, 1), dtype=uint8)
SCL_CLOUD_SHADOW_RES: numpy.ndarray(shape=(2, 4096, 4096, 1), dtype=uint8)
dataMask_RES: numpy.ndarray(shape=(2, 4096, 4096, 1), dtype=uint8)
}
scalar={
earth_sun_dist: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B02: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B03: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B04: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B05: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B06: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B07: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B08: numpy.ndarray(shape=(2, 1), dtype=float64)
}
)
# plot differences due to noise
tidx = 1
_, axes = plt.subplots(ncols=3, nrows=1, figsize=(15, 15), sharex=True, sharey=True)
axes[0].imshow(eop.data["BANDS-RAD-PAN_SHIFTED"][t_idx][..., [2, 1, 0]] / 100.0)
axes[0].set_title(f"BEFORE - {timestamp}")
axes[1].imshow(eop.data["L_out_RES"][t_idx][..., [2, 1, 0]] / 100.0)
axes[1].set_title(f"AFTER - {timestamp}")
axes[2].imshow(
np.abs(
eop.data["L_out_RES"][t_idx][..., [2, 1, 0]]
- eop.data["BANDS-RAD-PAN_SHIFTED"][t_idx][..., [2, 1, 0]]
)
)
axes[2].set_title("RGB difference")
plt.tight_layout()
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
8. Simulate Point Spread Function (PSF) Filtering¶
Step that simulates loss of signal due to limited bandwidth of the Module Transfer Function (MTF), which is equivalent to simulate the Point Spread Function (PSF) for the entire system (imager and platform).
The same executable as for SNR is used. If you have not downloaded it already, do so for your OS from this link. Once downloaded, allow the operating system to run the file and make it executable (e.g. chmod a+x phisat2_unix.bin). Then set the path to the executable in the cell below.
psf_filter_task = PhisatCalculationTask(
input_feature=(FeatureType.DATA, "L_out_RES"),
output_feature=(FeatureType.DATA, "L_out_PSF"),
executable=EXEC_PATH,
calculation="PSF",
)
eop = psf_filter_task(eop)
# check result:
eop
EOPatch(
bbox=BBox(((358232.0, 4679922.0), (377688.0, 4699378.0)), crs=CRS('32633'))
timestamps=[datetime.datetime(2021, 9, 2, 10, 8, 50), datetime.datetime(2021, 9, 7, 10, 8, 54)]
meta_info={
Shifts: {0: [(-0.1312835857123827, 0.6266645003616177), ...]<length=8>, ...}<length=2>
}
data={
BANDS-RAD-PAN_SHIFTED: numpy.ndarray(shape=(2, 4096, 4096, 8), dtype=float64)
L_out_PSF: numpy.ndarray(shape=(2, 4096, 4096, 8), dtype=float64)
L_out_RES: numpy.ndarray(shape=(2, 4096, 4096, 8), dtype=float64)
sunZenithAngles_RES: numpy.ndarray(shape=(2, 4096, 4096, 1), dtype=float32)
}
mask={
SCL_CIRRUS_RES: numpy.ndarray(shape=(2, 4096, 4096, 1), dtype=uint8)
SCL_CLOUD_RES: numpy.ndarray(shape=(2, 4096, 4096, 1), dtype=uint8)
SCL_CLOUD_SHADOW_RES: numpy.ndarray(shape=(2, 4096, 4096, 1), dtype=uint8)
dataMask_RES: numpy.ndarray(shape=(2, 4096, 4096, 1), dtype=uint8)
}
scalar={
earth_sun_dist: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B02: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B03: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B04: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B05: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B06: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B07: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B08: numpy.ndarray(shape=(2, 1), dtype=float64)
}
)
# visualize part of EOPatch before/after applying PSF kernel
n_timestamps = len(eop.timestamp)
fig, axes = plt.subplots(ncols=3, nrows=1, figsize=(15, 5))
axes[0].imshow(eop.data["L_out_RES"][0][2000:2200, 2000:2200, [2, 1, 0]] / 100)
axes[0].set_title(timestamp)
axes[1].imshow(eop.data["L_out_PSF"][0][2000:2200, 2000:2200, [2, 1, 0]] / 100)
axes[1].set_title("with PSF applied")
diff = eop.data["L_out_PSF"] - eop.data["L_out_RES"]
axes[2].hist(diff.ravel(), bins=40, log=True)
axes[2].set_title("difference, all bands and timestamps")
plt.tight_layout()
/var/folders/wv/d2zg_qp156z_4hsjy3yb2shm0000gn/T/ipykernel_95471/781997810.py:3: EODeprecationWarning: The attribute `timestamp` is deprecated, use `timestamps` instead. n_timestamps = len(eop.timestamp) Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
eop
EOPatch(
bbox=BBox(((358232.0, 4679922.0), (377688.0, 4699378.0)), crs=CRS('32633'))
timestamps=[datetime.datetime(2021, 9, 2, 10, 8, 50), datetime.datetime(2021, 9, 7, 10, 8, 54)]
meta_info={
Shifts: {0: [(-0.1312835857123827, 0.6266645003616177), ...]<length=8>, ...}<length=2>
}
data={
BANDS-RAD-PAN_SHIFTED: numpy.ndarray(shape=(2, 4096, 4096, 8), dtype=float64)
L_out_PSF: numpy.ndarray(shape=(2, 4096, 4096, 8), dtype=float64)
L_out_RES: numpy.ndarray(shape=(2, 4096, 4096, 8), dtype=float64)
sunZenithAngles_RES: numpy.ndarray(shape=(2, 4096, 4096, 1), dtype=float32)
}
mask={
SCL_CIRRUS_RES: numpy.ndarray(shape=(2, 4096, 4096, 1), dtype=uint8)
SCL_CLOUD_RES: numpy.ndarray(shape=(2, 4096, 4096, 1), dtype=uint8)
SCL_CLOUD_SHADOW_RES: numpy.ndarray(shape=(2, 4096, 4096, 1), dtype=uint8)
dataMask_RES: numpy.ndarray(shape=(2, 4096, 4096, 1), dtype=uint8)
}
scalar={
earth_sun_dist: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B02: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B03: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B04: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B05: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B06: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B07: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B08: numpy.ndarray(shape=(2, 1), dtype=float64)
}
)
# visualize part of EOPatch before/after applying PSF kernel
n_timestamps = len(eop.timestamp)
fig, axes = plt.subplots(ncols=3, nrows=1, figsize=(15, 5))
axes[0].imshow(eop.data["L_out_RES"][0][2000:2200, 2000:2200, [2, 1, 0]] / 100)
axes[0].set_title(timestamp)
axes[1].imshow(eop.data["L_out_PSF"][0][2000:2200, 2000:2200, [2, 1, 0]] / 100)
axes[1].set_title("with PSF applied")
diff = eop.data["L_out_PSF"] - eop.data["L_out_RES"]
axes[2].hist(diff.ravel(), bins=40, log=True)
axes[2].set_title("difference, all bands and timestamps")
plt.tight_layout()
/var/folders/wv/d2zg_qp156z_4hsjy3yb2shm0000gn/T/ipykernel_95471/781997810.py:3: EODeprecationWarning: The attribute `timestamp` is deprecated, use `timestamps` instead. n_timestamps = len(eop.timestamp) Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
9. Alternative Noise and PSF calculation¶
The AlternativePhisatCalculationTask allows users to input their own Signal To Noise ratios and their own Point Spread Function kernel values for addition of SNR and PSF.
kernel_bands = ["B1", "B2", "B3", "B0", "B7", "B4", "B5", "B6"]
psf_kernel = { band: np.random.random(size=(7,7)) for band in kernel_bands }
snr_bands = ["B02", "B03", "B04", "PAN", "B08", "B05", "B06", "B07"]
snr_values = { band: np.random.randint(20,250) for band in snr_bands }
snr_psf_task = AlternativePhisatCalculationTask(
input_feature=(FeatureType.DATA, "BANDS-RAD-PAN_SHIFTED"),
snr_feature=(FeatureType.DATA, "L_out_RES"),
psf_feature=(FeatureType.DATA, "L_out_PSF"),
snr_values=snr_values,
psf_kernel=psf_kernel,
l_ref=100
)
# uncomment to run
# eop = snr_psf_task(eop)
# check result:
eop
EOPatch(
bbox=BBox(((358232.0, 4679922.0), (377688.0, 4699378.0)), crs=CRS('32633'))
timestamps=[datetime.datetime(2021, 9, 2, 10, 8, 50), datetime.datetime(2021, 9, 7, 10, 8, 54)]
meta_info={
Shifts: {0: [(-0.1312835857123827, 0.6266645003616177), ...]<length=8>, ...}<length=2>
}
data={
BANDS-RAD-PAN_SHIFTED: numpy.ndarray(shape=(2, 4096, 4096, 8), dtype=float64)
L_out_PSF: numpy.ndarray(shape=(2, 4096, 4096, 8), dtype=float64)
L_out_RES: numpy.ndarray(shape=(2, 4096, 4096, 8), dtype=float64)
sunZenithAngles_RES: numpy.ndarray(shape=(2, 4096, 4096, 1), dtype=float32)
}
mask={
SCL_CIRRUS_RES: numpy.ndarray(shape=(2, 4096, 4096, 1), dtype=uint8)
SCL_CLOUD_RES: numpy.ndarray(shape=(2, 4096, 4096, 1), dtype=uint8)
SCL_CLOUD_SHADOW_RES: numpy.ndarray(shape=(2, 4096, 4096, 1), dtype=uint8)
dataMask_RES: numpy.ndarray(shape=(2, 4096, 4096, 1), dtype=uint8)
}
scalar={
earth_sun_dist: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B02: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B03: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B04: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B05: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B06: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B07: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B08: numpy.ndarray(shape=(2, 1), dtype=float64)
}
)
10. Compute L1C¶
If the specified processing level is L1C, reflectances are computed from radiances. Otherwise, radiances are provided as output for L1A and L1B.
reflectance_task = CalculateReflectanceTask(
input_feature=(FeatureType.DATA, "L_out_PSF"),
output_feature=(FeatureType.DATA, "PHISAT2-BANDS"),
processing_level=PROCESSING_LEVEL,
)
eop = reflectance_task(eop)
eop
EOPatch(
bbox=BBox(((358232.0, 4679922.0), (377688.0, 4699378.0)), crs=CRS('32633'))
timestamps=[datetime.datetime(2021, 9, 2, 10, 8, 50), datetime.datetime(2021, 9, 7, 10, 8, 54)]
meta_info={
Shifts: {0: [(-0.1312835857123827, 0.6266645003616177), ...]<length=8>, ...}<length=2>
}
data={
BANDS-RAD-PAN_SHIFTED: numpy.ndarray(shape=(2, 4096, 4096, 8), dtype=float64)
L_out_PSF: numpy.ndarray(shape=(2, 4096, 4096, 8), dtype=float64)
L_out_RES: numpy.ndarray(shape=(2, 4096, 4096, 8), dtype=float64)
PHISAT2-BANDS: numpy.ndarray(shape=(2, 4096, 4096, 8), dtype=float64)
sunZenithAngles_RES: numpy.ndarray(shape=(2, 4096, 4096, 1), dtype=float32)
}
mask={
SCL_CIRRUS_RES: numpy.ndarray(shape=(2, 4096, 4096, 1), dtype=uint8)
SCL_CLOUD_RES: numpy.ndarray(shape=(2, 4096, 4096, 1), dtype=uint8)
SCL_CLOUD_SHADOW_RES: numpy.ndarray(shape=(2, 4096, 4096, 1), dtype=uint8)
dataMask_RES: numpy.ndarray(shape=(2, 4096, 4096, 1), dtype=uint8)
}
scalar={
earth_sun_dist: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B02: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B03: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B04: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B05: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B06: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B07: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_B08: numpy.ndarray(shape=(2, 1), dtype=float64)
sol_irr_PAN: numpy.ndarray(shape=(2, 1), dtype=float64)
}
)
# plot difference dub to radiance to refelectance conversion if applied
t_idx = 1
viz_factor = 3.5 if PROCESSING_LEVEL.value == ProcessingLevels.L1C.value else 0.01
_, axes = plt.subplots(ncols=2, nrows=1, figsize=(15, 15), sharex=True, sharey=True)
axes[0].imshow(eop.data["L_out_PSF"][t_idx][..., [2, 1, 0]] / 100.0)
axes[0].set_title(f"L_out_PSF - {timestamp}")
axes[1].imshow(eop.data["PHISAT2-BANDS"][t_idx][..., [2, 1, 0]] * viz_factor)
axes[1].set_title(f"PHISAT2-BANDS - {timestamp}")
plt.tight_layout()
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
remove_feature_task3 = RemoveFeatureTask(
[
(FeatureType.DATA, "BANDS-RAD-PAN_SHIFTED"),
(FeatureType.DATA, "L_out_PSF"),
(FeatureType.DATA, "L_out_RES"),
]
)
eop = remove_feature_task3(eop)
11. Gridding¶
On board of Φ-sat-2 the images covering the swath width will be split into smaller patches for computational processing. For creation of the training dataset, non-overlapping image patches of size 256x256 are created. When testing your algorithm, the image patches should possibly overlap to remove edge artefacts.
Change the parameters in the cell below to modify the size of the image patches and their overlap.
CELL_SIZE = 256
GRID_OVERLAP = 0.0
The gridding task works for a single timestamp at a time. If you needed to grid multiple timeframes, you need to have a task for each timeframe to grid and export to Tiff.
time_index = 0
FEATURE_NAMES = [
(FeatureType.DATA, "PHISAT2-BANDS"),
(FeatureType.MASK, "SCL_CIRRUS_RES"),
(FeatureType.MASK, "SCL_CLOUD_RES"),
(FeatureType.MASK, "SCL_CLOUD_SHADOW_RES"),
(FeatureType.MASK, "dataMask_RES"),
]
gridding_tasks = []
for feature_type, feature_name in FEATURE_NAMES:
gridding_tasks.append(
GriddingTask(
(feature_type, feature_name),
(feature_type, f"{feature_name}-GRID_{time_index}"),
(FeatureType.VECTOR_TIMELESS, f"{feature_name}-GRID_{time_index}"),
size=CELL_SIZE,
overlap=GRID_OVERLAP,
resolution=PHISAT2_RESOLUTION,
time_index=time_index,
)
)
eop = gridding_tasks[0](eop)
NSAMPLES = 9
_, axs = plt.subplots(
figsize=(15, 15), ncols=3, nrows=NSAMPLES // 3, sharex=True, sharey=True
)
for nsample in range(NSAMPLES):
axs[nsample % 3][nsample // 3].imshow(
eop.data[f"PHISAT2-BANDS-GRID_{time_index}"][nsample][..., [2, 1, 0]]
* viz_factor
)
plt.tight_layout()
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
# show resulting grid
fig, axs = plt.subplots(figsize=(15, 15))
axs.imshow(
eop.data["PHISAT2-BANDS"][1][..., [2, 1, 0]] * viz_factor, extent=get_extent(eop)
)
eop.vector_timeless[f"PHISAT2-BANDS-GRID_{time_index}"].boundary.plot(ax=axs);
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
12. Export to Tiff¶
Creates tiffs for a single timestamp. The simulated bands along with the cloud, cloud shadow, and cirrus mask are exported as tiff files. These files can be visualised and further analysed into any GIS software, such as QGis.
The format of the tiff filename is as follows:
{utm_easting_centroid}-{utm_northing_centroid}_{epsg_crs}_{feature}_{datetime}_{cell_number}
# create a folder for output tiffs
os.makedirs("tiff_folder", exist_ok=True)
export_tiff_tasks = []
for feature_type, feature_name in FEATURE_NAMES:
export_tiff_tasks.append(
ExportGridToTiff(
data_stack_feature=(feature_type, f"{feature_name}-GRID_{time_index}"),
grid_feature=(
FeatureType.VECTOR_TIMELESS,
f"{feature_name}-GRID_{time_index}",
),
out_folder="./tiff_folder",
time_index=time_index,
)
)
eop = export_tiff_tasks[0](eop)
# check result:
os.listdir("tiff_folder")[:10]
['367960-4689650_32633_PHISAT2-BANDS-GRID_0_2021-09-02T10-08-50_000.tiff', '367960-4689650_32633_PHISAT2-BANDS-GRID_0_2021-09-02T10-08-50_001.tiff', '367960-4689650_32633_PHISAT2-BANDS-GRID_0_2021-09-02T10-08-50_002.tiff', '367960-4689650_32633_PHISAT2-BANDS-GRID_0_2021-09-02T10-08-50_003.tiff', '367960-4689650_32633_PHISAT2-BANDS-GRID_0_2021-09-02T10-08-50_004.tiff', '367960-4689650_32633_PHISAT2-BANDS-GRID_0_2021-09-02T10-08-50_005.tiff', '367960-4689650_32633_PHISAT2-BANDS-GRID_0_2021-09-02T10-08-50_006.tiff', '367960-4689650_32633_PHISAT2-BANDS-GRID_0_2021-09-02T10-08-50_007.tiff', '367960-4689650_32633_PHISAT2-BANDS-GRID_0_2021-09-02T10-08-50_008.tiff', '367960-4689650_32633_PHISAT2-BANDS-GRID_0_2021-09-02T10-08-50_009.tiff']
nodes = linearly_connect_tasks(
scl_download_task,
scl_cloud_task,
input_task,
# filter_task, # Uncomment this if you want to filter by 100% of valid data
add_meta_task,
radiance_task,
add_pan_task,
casting_task,
*resize_task_list,
remove_feature_task1,
band_shift_task,
remove_feature_task2,
crop_task,
snr_task,
psf_filter_task,
reflectance_task,
remove_feature_task3,
*gridding_tasks,
*export_tiff_tasks,
)
workflow = EOWorkflow(nodes)
# This line repeats the simulation as done in the cells above
# results = workflow.execute({nodes[0]: {"bbox": bbox, "time_interval": time_interval}})
14. Core Dataset¶
To get you started with your AI application at-the-edge, we simulated Φ-sat-2 imagery for around 490 location on Earth. The workflow as defined above was used to simulate L1C products for the locations and dates as specified in the provided phisat2-locations.gpkg file. The code below was used to generate such images.
The locations and acquisition times seek to match the ones for the PRISMA images utilized in the IMAGIN-e track of the challenge. The information about the time difference from the PRISMA acquisition time is stored in the GeoDataFrame, as well as the maxcc of Sentinel-2 image and the overlap between the Φ-sat-2 bounding box and the geometry of the Sentinel-2 image. The dataframe also contains the name of the coresponding PRISMA image prisma_name, in case you wanted to cross-reference the datasets.
We suggest you to prepare a similar dataframe with dates and locations suitable for your application and use the code below to simulate Φ-sat-2 imagery.
For the locations provided we suggest you download directly the provided tiff files, unless you want to tweak the simulation workflow. The name of the zip archive containing the image patches can be found in the zipfile column, so you can retrieve from storage only the simulated locations of interest. The zipfiles are available on cloudferro on the following public bucket
https://s3.waw2-1.cloudferro.com/swift/v1/AUTH_afccea586afd4ef3bb11fe37dd1ddfbf/OrbitalAI_Data/,
and you can download any file by appending hte zip archive name, e.g.
phisat2_locations = gpd.read_file("phisat2-locations.gpkg")
len(phisat2_locations)
489
fig, ax = plt.subplots(figsize=(15, 10))
WORLD_GDF.plot(ax=ax, color="b", alpha=0.2)
WORLD_GDF.boundary.plot(ax=ax, color="b")
phisat2_locations.plot(ax=ax, color="r");
phisat2_locations.head(5)
| datetime_s2 | intersection_ratio | cloud_cover_s2 | time_difference | prisma_name | zipfile | geometry | |
|---|---|---|---|---|---|---|---|
| 0 | 2019-10-12T16:28:28 | 1.000000 | 0.30 | -14 | PRS_L1_STD_OFFL_20191026162600_20191026162604_... | 335040-5955400_32618_2019-10-12T16-28-28.zip | POINT (-77.50000 53.72120) |
| 1 | 2019-11-30T15:30:28 | 1.000000 | 16.20 | -7 | PRS_L1_STD_OFFL_20191207152456_20191207152500_... | 632700-1080910_32618_2019-11-30T15-30-13.zip | POINT (-73.79000 9.77629) |
| 2 | 2019-12-18T08:48:42 | 1.000000 | 28.67 | -1 | PRS_L1_STD_OFFL_20191219085427_20191219085431_... | 539510-4571900_32636_2019-12-18T08-48-42.zip | POINT (33.47200 41.29760) |
| 3 | 2019-12-20T09:31:05 | 0.805958 | 3.45 | -1 | PRS_L1_STD_OFFL_20191221093011_20191221093015_... | 495240-3696340_32634_2019-12-20T09-31-01.zip | POINT (20.94890 33.40640) |
| 4 | 2019-12-26T09:50:55 | 1.000000 | 0.23 | 4 | PRS_L1_STD_OFFL_20191222094637_20191222094641_... | 688270-3825460_32633_2019-12-26T09-50-55.zip | POINT (17.05200 34.55370) |
fig, ax = plt.subplots(figsize=(20, 5), ncols=3)
phisat2_locations.intersection_ratio.hist(ax=ax[0], bins=20, log=True)
ax[0].set_title("Overlap between Φ-sat-2 bbox and S2 image")
ax[0].set_xlabel("Overlap fraction")
ax[0].set_ylabel("log (# locations)")
phisat2_locations.time_difference.hist(ax=ax[1], bins=20)
ax[1].set_title("Time difference between S2 and PRISMA")
ax[1].set_xlabel("# days")
ax[1].set_ylabel("# locations")
phisat2_locations.cloud_cover_s2.hist(ax=ax[2], bins=20)
ax[2].set_title("Cloud cover of S2 images")
ax[2].set_xlabel("maxcc")
ax[2].set_ylabel("# locations");
The simulation workflow can be parallelized over the locations to speed up the processing. EOExecutor allows to scale the processing over multiple processors/instances. For this to work we need to prepare the list of arguments for the executor, and the number of processing units to parallelize the execution over.
exec_args = [
{
nodes[0]: {
"bbox": get_utm_bbox(coords.y, coords.x),
"time_interval": (atime.split("T")[0], atime.split("T")[0]),
}
}
for atime, coords in zip(phisat2_locations.datetime_s2, phisat2_locations.geometry)
]
exec_args[0]
{EONode(task=<eolearn.io.sentinelhub_process.SentinelHubInputTask object at 0x7f07cde365c0>, inputs=(), name='SentinelHubInputTask'): {'bbox': BBox(((324970.0, 5945320.0), (345110.0, 5965460.0)), crs=CRS('32618')),
'time_interval': ('2019-10-12', '2019-10-12')}}
eoexecutor = EOExecutor(workflow=workflow, execution_kwargs=exec_args, save_logs=True)
eoexecutor.run(workers=4);